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LÉVY LAWS FOR LATTICE GAS AUTOMATA 
OLIVIER TRIBEL AND JEAN PIERRE BOON 


ABSTRACT. We consider the use of power law distributions P(z) ~| z |-17^ (1 < 
и < 2) for lattice gas automata (LGA). 

(а) We implement the distribution P(£) for particle displacements with length 
|£ | in a LGA and we investigate the resulting diffusive transport. We show 
that anomalous diffusion (var(| £ |) ~ 12/1} cannot be observed, but that normal 
diffusion ((/?) ~ Dt) can be tuned so that the diffusion coefficient D takes values 
larger than those obtained from usual random walk procedures. 

(b) We present а new LGA model with interaction range r set according to 
P(r); we show that the static properties of the LGA depend on Efr] whereas the 
dynamic properties depend on higher moments of the distribution. 
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1. REMARKS ON STABLE LAWS AND DIFFUSION. 
Consider the power law distribution 
Pla) ^| z [UU Ve 2 0, (1) 
which belongs to the basin of attraction of the stable (Lévy) law L,g(x), where 
B measures the asymmetry of the distribution (here 3 = 0) and u determines the 


behavior of the distribution for large values of z. For u < 2, the moments are such 
that 


(2) 


Де") forv > n. 


tae for v < и, 


The normal distribution corresponds to и = 2, 8 = 0 and any distribution defined 
on a bounded domain belongs to the basin of attraction of the normal distribution. 

If S, = X, Xi has the same distribution as the equally distributed random 
variables X;, the probability distribution P(S,,) is said to be stable. АП probability 
densities of laws that belong to the basin of attraction of a given stable law share the 
same behavior at infinity, which is identical to the behavior of the probability density 
of the stable law itself. However, practical implementations are not performed with 
distribution functions, but with realizations of random sampling, i.e. in practice, 
the limit n — oo is unrealistic. What are the consequences of the limitation of the 
sampling domain? 

(i) Consider the sum Sy = Y. Xj, where the N terms are distributed according 
to the power law (1). Obviously, the theoretical sampling domain (N — oo) can- 
not be probed completely. However what can be determined is the most probable 
largest z defined such that a larger value will be obtained only once amongst the N 
realizations, 1.e. 


P(z;) = 一 3 
PELLE. (3) 
wherefrom it follows that £, ~ МИР. Consequently, Sw is set by its largest term, 
and 


E[Sy] = NELX] = N s 27 о NUS (4) 


As a result, there exists a natural cutoff le which bounds the distribution. Note that 
for 0 < < 1, le grows faster than N, which precludes the existence of any moment 
of Sw. For 1 < и < 2, and for sufficiently large N, the sampling procedure can be 
performed, and E[Sy] exists. 

(ii) If (independently of the criteria for stable laws) an upper bound £,, is imposed 
to the distribution density, then the distribution belongs to the basin of attraction 
of the normal law. Hence, if N is sufficiently large, Sx is normally distributed and 
E[Sy] and var[Sw] exist. In conclusion, with an imposed cutoff ба (< (.) and a given 
value of µε [1,2], а normal distribution can be generated with a power law. 
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(iii) The characteristic function of the distribution is defined as 


E(k) = Ее | = > e7 P(x) (5) 


Consider a random walk generated by successive jumps of variable length at dis- 
crete time steps. After N jumps (or equivalently after a time t = N time steps), 
the position of the random walker is 


N 
i—1 


where L; is the length of the jump at the ¿th step. The (в are distributed according 
to the power law Р(6) which belongs to the domain of attraction of the stable (Lévy) 
law L olx), the characteristic function of which is 


E(k) = е7", (7) 
(a) u = 2 corresponds to the normal distribution with variance 
ха X] = N var|/] = N2dD = 2d Dt, (8) 


where d is the space dimension and D is the diffusion coefficient. (8) is the expected 
result for normal diffusion. 

(b) for 1 < u < 2, and with a natural cutoff le ~ N'/" = H/H (see above), one 
finds 


£c 
var[X] ~ Sra? Lio(x) ~ ple, (9) 
æ=1 


This result expresses that the mean square displacement grows faster than t, a 
characteristic feature of anomalous diffusion (see Bouchaud and George|1990]). 


2. LATTICE GAS DIFFUSION 


Consider a lattice gas automaton (LGA) where particles diffuse on a lattice £ 
through the iterated action of two successive operators: 
(i) The propagation operator P displaces particles from one node to another along 
the direction of the particle's velocity vector ci (¿ = 1,... „6 denotes the channel 
occupied by the particle, and there are b channels with velocity modulus | e; | — 1). 
(ii) The particles residing on a node are redistributed amongst the channels through 
the action of an operator R. The sequence R o P takes one discrete time step. The 
Boolean variable n?(r,t) denotes the occupation of channel ? at the node located at 
position r(€ £) at time t. The propagation step is formally expressed by 


ni(r,t):PW(r— (ei, t£ — 1) = $ C ni(r — (ei t 1) (10) 
2 


where Ç is a Boolean variable with power law distribution P(4), and í (=integer) is 
the length of the displacement along the direction ci. The sum on the RHS of (10) 
has a lower bound £ = 1; the upper bound will be discussed below. 
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The redistribution of channel occupations is operated by 


n; i, t): RAS r, t) = Xi nil r, t); (1-4 j. mod b), (11) 


where 6; is a Boolean variable equally distributed over the b channels. 
Defining the occupation probability as 


n(r,e t) =< ni(r,t) >, (12) 
one obtains from Eqs. (10) and ul the lattice Boltzmann equation which reads 


n(r,c;; t4 1) = EEPO (r — lej, cj; t). (13) 


Introducing the space Fourier transform of n(r, ci; t): 
nelci t) = nkelkr tot, (14) 


we obtain from (13) and (14) the eigenvalue equation: 
。 1 
Уна) (15) 
7 


with k; = k - cj, and where £(k) is the characteristic function of the distribution of 
jumps (see (7)). 

We are interested in the long wavelength limit, i.e. klm < 1. Two cases arise: 
when the distribution of jumps belongs to the basin of attraction of a stable law, i.e. 
the summation in (10) extends to infinity, the dispersion equation obtained from 


(15) and (7) reads 
wlk) > ΙΡ. (16) 


When the upper bound of the summation in (10) is set to a finite value lm, or when 
the power parameter и = 2, one has the normal diffusion result, namely 


2 2 
w(k) = ビー —Dk (17) 
where с? is the variance of the distribution (see (8)). 

The usual LGA implementation for diffusion proceeds with single lattice unit 
length displacements, which in the present description restricts the summation in 
(10) to the term £ = 1. We then retrieve the usual result D = По = (2d)! (= 1/4 
in two-dimensional space) (Dab εἰ al.[1993]). 

We performed simulations on a square lattice with displacements distributed ac- 
cording to the power law P(£) ~ СГ for 1 € £ € ἐν (and P(£) = 0 otherwise). We 
measured the mean square displacement as a function of the number of time steps; 
typical results are shown in Fig. 1. Besides the expected result for normal diffusion 
(Do = 0.25, when £ = 1), two interesting features are obtained: (i) provided lm 
is not too large, normal diffusion holds (i.e. (2) is a linear function of t), but the 
diffusion coefficient increases (D > Do) with lm and with и (1 < u < 2) (a typical 
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illustration is given in Fig. 1.a); (ii) for large values of lm (typically 4m = O(L), 
where L is the linear size of the lattice), the predicted 1?/" behavior cannot be ob- 
served (see Fig. 1.b). Another set of simulations was conducted under ”source-sink” 
conditions: particles are injected at one boundary of the lattice and annihilated at 
the opposite boundary. The corresponding solution to the diffusion equation for the 
particle density reads 


r 

n(r,t) = по е с (==) (18) 
which asymptotically yields a steady state solution linear in r. А typical result is 
shown in Fig. 2.a. (lm = 5, и = 1.5). Excellent agreement is obtained between the 
simulation data and the theoretical prediction. In the absence of imposed cutoff, 
the macroscopic equation governing the diffusion process is unknown; the simulation 
results shown in Fig. 2.b (и = 1.5 and ба, = 256 = L) cannot be computed by the 
theory. This result, as well as the preclusion to observe anomalous diffusion (+ {24} 
in a bounded CA universe, confirms that real long tailed distributions cannot be 
implemented in finite systems using integers. 


3. LATTICE GAS AUTOMATA WITH INTERACTION RANGE 


We present a LGA with non-local interactions which is a generalization of the 
model constructed by Appert, Zaleski, and d'Humieres [1993], and further analyzed 
by Gerits, Ernst, and Frenkel [1993]. The basic geometry is the model with hexag- 
onal symmetry initially proposed by Frisch, Hasslacher, and Pomeau [1986]. Аз 
illustrated Im Fig. 3, particles can interact non locally by exchanging momentum at 
a distance r, r being distributed according to a power law P(r). The automaton 
evolution is governed by an iterated sequence of three operators: PoCoZ, where P is 
the propagation operator (see Eq.(10), with only one term in the sum), C is the local 
collision operator (see e.g. Frisch et al.[1987]), and Z is the non-local interaction 
operator 


210)=5 У 


Here, the channels involved are those with | с |= 1, (1 = 0,...,5) and (¿i + 7) is 
taken modulo 6, whereas the local collision operator redistributes particles amongst 


pv) У i] (19) 


|= 


all channels including the zero-velocity channel | св |= 0. /'(п*) describes explicitly 
the non-local interactions, e.g. for the configuration change (asco) — (aocs) of Fig.3, 
one has 


K (n) = TG A ni). ^ (οὔ Λ ni) "i V [c ^ nj). ^ (ni ^ 3) (20) 


In the long wavelength limit and for homogeneous density, one can derive the com- 
pressibility equation (see e.g. Gerits εἰ al.[1993]). From the Boltzmann equation, 


6 OLIVIER TRIBEL AND JEAN PIERRE BOON 


one obtains 


9 2 cin(x,ci;t) = —V : P. (21) 


Р = M cici n(x, cj; t ;X Y ССр} 2 rP(r) Ir i (n). (22) 
i i j=-1 

On the RHS of (22), the first term is the kinetic contribution to the pressure tensor 

P and the second term arises from non-local interactions. At global equilibrium, 

the distribution function is n = p/7, where p is the average density per node, and 

it follows from (22) that the hydrostatic pressure is given by 


p= 5e (1-0) Xr Per (23) 


where we observe that the interactive part depends on the first moment ( of the 
distribution. The compressibility equation (23) is illustrated in Fig. 4 for various 
values of (г). From Eq.(23), one observes that the compressibility p^! (Op/Op) or the 
square of the sound velocity c? = др/др becomes negative when (г) exceeds some 
critical value r.(~ 8). As a result pressure fluctuations would critically amplify: 
acoustic mode explosion is prevented by a regime transition as the system undergoes 
phase separation when the average interaction range exceeds the critical value re (see 
Fig. 5). The evolution of the density distribution during phase separation is shown 
in Fig. 6. 

We now consider the dynamic properties of the interactive LGA. We essentially 
follow the standard procedure (see Grosfils εἰ al.[1993], Gerits et al.[1993]) which 
is applied here to the probabilistic dynamics of the LGA with interaction range 
governed by the distribution P(r). We linearize the Boltzmann equation to obtain 


ón(x + ει. οι; Ὁ) = (1 + Q)ia ón (x, ca; t), (24) 


ón'(x,cq;t) = (1 + A(r))qrdn(x, си; t) (25) 


where Q; is the linearized collision operator and A(r),¢ is the linearized interaction 
operator. By eigenmode decomposition 


n(x, ci; t) = V, (К, ci) exp (ik x + za(k)t) (26) 
we obtain the eigenvalue equation 
[eer — (1 + NA + ACK) V, (k,c4) = 0 (27) 
with 
A (k)ón, = τΣθω er F(\)G(6\), (28) 


where the functions Z and G involve products of n (the global equilibrium distri- 
bution function n = p/7) and of ón's taken at interaction distances. Schematically, 
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the result can be written in the form 


A; (k)ón, ~ — ralni 一 6n4)(1 — Е(К)) + ска(дп + 67。) (k) (29) 
with the static correlation functions 
ко = п(1 — n), кз = (1 — 2n)&a (30) 


and where £&'(k) and £”(k) are the real and imaginary parts of the characteristic 
function. 

Consequently it follows that the dynamic properties, in particular the dynamic 
structure factor, which can be computed from the linearized Boltzmann equation, 
will depend on the distribution governing the microscopic dynamics. The transport 
coefficients can be obtained by considering the eigenvalue problem for the hydrody- 
namic modes in the long wavelength limit, i.e. by low Е expansion. This requires 
that the first moments of the distribution exist and that, provided 4E[r], 3 var[r],..., 
k(r), k?(r?),... be small (< 1). The details of the analysis are presented elsewhere 
(Tribel and Boon[1993]). 

In conclusion, we have shown that when probabilistic dynamics is used to model 
LGA's with long range interactions, the static properties involve only the first mo- 
ment of the distribution, and, in that respect, no essential difference is found with 
minimal models with fixed interaction distance; however, we give indication that 
the dynamic properties will depend on higher moments of the distribution, a feature 
which should meet the well-known sensitivity of the transport coefficients to the 
interaction potential in real systems. 
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Figure Captions. 


Fig. 1. Mean square displacement as a function of time. (а) Normal diffusion: 
simulation data (black dots) and linear fit (solid lines). Usual normal diffusion: 
(= 1, Do = 0.25 (lower set of data), and enhanced normal diffusion: m = 5, и = 
1.8, D = 0.5 (upper set of data). (b) "Abnormal" diffusion: £,, = 512, и = 1.5, 
five sets of data with 1024 test particles. The mean square displacement is given in 
lattice unit lengths and the time in number of automaton time steps. Lattice size: 
512*512; periodic boundary conditions. 

Fig. 2. Density profile distribution in "source-sink" simulation. (a) lm = 5, и = 1.5: 
density profile after 6580 automaton time steps. Simulation data (noisy curve) and 
theoretical curve, Eq. (18). (b) а = 256, = 1.5: simulation results for density 
profiles at successive times up to 10? automaton time steps. The density is expressed 
in particles per channel as a function of node position |0, 255]. 

Fig. 3. Attraction interactions on a triangular lattice. Arrows (double arrows) 
represent channel occupations before (after) interaction between nodes separated 
by r lattice unit lengths. Configurational changes (або) — (aob1), (asco) 一 
(аосз), (asdo) 一 (aod5) occur with equal probability (= 1/3). 

Fig. 4. Pressure as a function of density, Eq. (23), for various values of (r) (—5, 8, 
12, 20 from upper to lower curve). The density n is given in particle number per 
channel. 

Fig. 5. Successive states of the automaton during phase separation. Time intervals 
during states are 100-200 time steps. r is distributed uniformly over the range [1,20] 
with (r) = 10.5. The initial state is homogeneous. Black (white) dots represent 
nodes with density larger (smaller) than the average density (n — 0.25). Lattice 
size: 256*256. 

Fig. 6. Same as Fig.5 showing the evolution of the density distribution during phase 
separation, starting from initial homogeneous state (n — 0.25). 
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